C*********************************************************************C
C* *C
C* olchem.for (olchem.msf on mplvax.sri.com) *C
C* *C
C* Written by: David L. Huestis, Molecular Physics Laboratory *C
C* *C
C* Copyright (c) 1981-84,1989,1991,1998 SRI International *C
C* All Rights Reserved *C
C* *C
C* This software is provided on an as is basis; without any *C
C* warranty; without the implied warranty of merchantability or *C
C* fitness for a particular purpose. *C
C* *C
C*********************************************************************C
C
C EDIT HISTORY:
C
C 05-23-98 DLH Filename/portability repairs and documentation
C
C 09-04-91 DLH First attempt at portable version, not used
C
C 04-04-89 DLH Increase number of species to 60
C
C 03-24-84 DLH Last unrecorded changes on PDP-11/40
C
C 02-08-84 DLH MAKE A FILE OF REACTION RATES CALLED OLCHEM.RXR
C
C 01-07-83 DLH RE-DIMENSION TO 200 REACTIONS, 40 SPECIES
C
C 01-23-82 DCL Don Lorents changed to write time and species
C concentrations in a file called OLCHEM.CON
C
C 11-23-81 DLH CONVERTED FOR MPL PDP 11/40 from RT-11 version
C obtained from Dave Golden.
C
C --------------------------------------------------------------------
C
PROGRAM OLCHEM ! program name of SRI International version
C PROGRAM CHEMK ! original program name
C
C ---- What follows is the documentation from the 1970s original ----
C
C CHEMK IS A FORTRAN PROGRAM WHICH, WHEN GIVEN A PREDETERMINED
C KINETIC MECHANISM, COMPUTES THE CONCENTRATIONS OF THE VARIOUS
C REACTANTS IN TIME. IT WAS WRITTEN BY GARY Z. WHITTEN OF SYSTEMS
C APPLICATIONS INCORPORATED IN SAN RAFAEL, CALIFORNIA AND INCORPORATES
C THE GEAR INTEGRATION PACKAGE OF A. HINDMARSH OF LAWRENCE LIVERMORE
C LABORATORIES IN LIVERMORE, CALIFORNIA. THE PRINTER-PLOT ROUTINE WAS
C PROVIDED BY DAVID C. WHITNEY OF SAI AND THE PROGRAM WAS CONVERTED TO
C ANSI FORTRAN BY JIM MEYER OF SAI. THE COMPUTER CODE THAT FOLLOWS
C IS EXECUTABLE ON ANY IBM OR CDC MACHINE WITH A FORTRAN IV COMPILER.
C
C TWO SEPARATE SETS OF DATA ARE NEEDED TO EXECUTE THE PROGRAM.ONE SET OF
C DATA CONTROLS THE COMPUTATIONAL PART OF THE PROGRAM AND IS
C NECESSARY TO ESTABLISH THE INTEGRATION ROUTINE. THE SECOND SET OF DATA
C CONTROLS THE PLOTTER OUTPUT AND PROVIDES THE PARAMETERS NECESSARY TO
C SPECIFY THE OUTPUT FORMAT. THE READER IS REFERED TO THE PROGRAM WRITE
C UP FOR FURTHER INFORMATION.
C
C GARY Z. WHITTEN
C SYSTEMS APPLICATIONS INCORPORATED
C 950 NORTHGATE DRIVE
C SAN RAFAEL, CALIFORNIA 94903
C (415) 472-4011
C
C --------------------------------------------------------------------
C
PARAMETER (MAXSPC=60, MSMS1=MAXSPC*(MAXSPC+1))
REAL IRS,IBLANK,MBLANK,JINTER,ITITLE,IGO
COMMON /DATA/ NR, KR( 200,7), A( 200), S( 200), ITITLE(7),
C TEMP, ERR,
C START, STOPP, PC, NP,
C SIG(60),IP(60), ITYPE(200), R( 200),
C BK, SG, DILUT
COMMON /NAMES/ SPECIS(60), REACT(60), NS
COMMON /STCOM1/ N, T, GUESS,
C HMIN, HMAX, EPS1,
C MF1, KFLAG1, JSTART
COMMON /STCOM2/ YMAX(60)
COMMON /STCOM3/ ERROR(60)
COMMON /STCOM4/ PW(MSMS1)
COMMON /STCOM5/ FSAVE( 120)
COMMON /ALPHA/ IBLANK,MBLANK,JINTER
COMMON /INOUT/ IN, IOUT, ITAPE,LUNRXR
DIMENSION C(360),IRS(200,7),IGO(2)
1 FORMAT(7A4)
2 FORMAT(14H THE REACTIONS)
3 FORMAT(I3,7A4,2E12.5)
4 FORMAT(/,I5,3X,3(1X,A4,1X),1H=,4(1X,A4,1X),1P2E12.4)
5 FORMAT(30X,7A4,//,23H INITIAL CONCENTRATION ,//
C (10X,5(4X,A4,4X)))
6 FORMAT(/(8X,1P 5E12.3))
7 FORMAT(I3,I6,3F10.0)
8 FORMAT(/,34H THE TEMPERATURE OF THE CELL WAS =,1PE9.2,
C 26H AND THE ERROR TOLERANCE =,E9.2)
9 FORMAT(/,29H THE RATE CONSTANTS USED WERE,/,(/,8X,1P 5E12.3))
10 FORMAT(/,30H THE OVERALL DILUTION RATE WAS,1PE9.2)
12 FORMAT(8F10.0)
13 FORMAT(3A4,48X,I5)
14 FORMAT(I3,A4,F10.0)
15 FORMAT(7E10.3)
C
C ALPHAMERIC DATA
C
DATA IBLANK/4H /,MBLANK/4HM /,JINTER/4HINTV/
DATA IGO/4HMORE,4HCONT/
C
C INPUT/OUTPUT DEVICES
C
IN=1
IOUT=5
ITAPE=2
LUNRXR=3
C
C$IF DEFINED( pdp11 )
C CALL ASSIGN(IOUT,'OLCHEM.LST')
C CALL FDBSET(IOUT,'UNKNOWN')
C CALL ASSIGN(IN,'OLCHEM.DAT')
C CALL FDBSET(IN,'READONLY')
C CALL ASSIGN(ITAPE,'OLCHEM.CON')
C CALL FDBSET(ITAPE,'UNKNOWN')
C CALL ASSIGN(LUNRXR,'OLCHEM.RXR')
C CALL FDBSET(LUNRXR,'UNKNOWN')
C
C$ELSEIF DEFINED( vms )
C OPEN(UNIT=IOUT,FILE='OLCHEM.LST',STATUS='UNKNOWN',
C 1 carriagecontrol='list')
C OPEN(UNIT=IN,FILE='OLCHEM.DAT',STATUS='OLD',READONLY)
C OPEN(UNIT=ITAPE,FILE='OLCHEM.CON',STATUS='UNKNOWN',
C 1 carriagecontrol='list')
C OPEN(UNIT=LUNRXR,FILE='OLCHEM.RXR',STATUS='UNKNOWN',
C 1 carriagecontrol='list')
C
C$ELSEIF DEFINED( unix )
C this is really for anything else (i.e. any FORTRAN-77)
OPEN(UNIT=IOUT,FILE='OLCHEM.LST',STATUS='UNKNOWN')
OPEN(UNIT=IN,FILE='OLCHEM.DAT',STATUS='OLD')
OPEN(UNIT=ITAPE,FILE='OLCHEM.CON',STATUS='UNKNOWN')
OPEN(UNIT=LUNRXR,FILE='OLCHEM.RXR',STATUS='UNKNOWN')
C
C$ELSE
C we assume that this is MS-FORTRAN
C OPEN(UNIT=IOUT,FILE='OLCHEM.LST',STATUS='UNKNOWN')
C OPEN(UNIT=IN,FILE='OLCHEM.DAT',STATUS='OLD',MODE='READ')
C OPEN(UNIT=ITAPE,FILE='OLCHEM.CON',STATUS='UNKNOWN')
C OPEN(UNIT=LUNRXR,FILE='OLCHEM.RXR',STATUS='UNKNOWN')
C$ENDIF
C
C CLEAR PATTERN MATRIX AND SET THE FIRST ELEMENTS
C
DO 20 J=1, 200
DO 20 K=1,7
20 KR(J,K)=0
DO 25 J=1, 200
A(J)=0.
R(J)=0.
25 S(J)=0.
NR=0
NS=1
80 NS=NS-1
C
C NX = LAST NUMBER OF THE REACTION SEQUENCE
C NPLOT = PLOT OPTION
C DILUT= DILUTION FACTOR
C
READ(IN,14)NX,NPLOT,DILUT
call do_ff(iout)
WRITE(IOUT,2)
C
C REACTION INPUT DATA
C
11 READ(IN,3)J,(IRS(J,I),I=1,7),A(J),S(J)
WRITE(IOUT,4)J,(IRS(J,K),K=1,7),A(J),S(J)
KR(J,1)=100
IF(J.GT.NR) NR=J
IF(J.NE.NX) GO TO 11
C
C ESTABLISH REACTION MATRIX
C
CALL MATRX (C,IRS)
C
C TITLE CARD AND OPTIONS
C
90 READ(IN,1)(ITITLE(I),I=1,7)
NFL=0
IF(ITITLE(1).EQ.IGO(1))GO TO 80
IF(ITITLE(1).EQ.IGO(2))GO TO 98
IF(ITITLE(1).EQ.IBLANK) STOP
C
35 READ(IN,7)N,NPRNT,TPRNT,TSTEP,TFACT
C
C TIME STEP SKIP OPTION
C
IF(NPRNT.EQ.0) NPRNT=10
C
C TIME STEP LENGTH OPTION
C
IF(TPRNT.EQ.0.)TPRNT=1.E+20
C
C CONCENTRATION OF SPECIES INITIALLY PRESENT
C
READ(IN,1)(REACT(I),I=1,N)
READ(IN,15)(C(I),I=1,N)
C
C STARTING AND ENDING INTEGRATION TIMES
C
READ(IN,12)START,STOPP,TEMP,ERR
C
C SPECIFICATION OF THE TEMPERATURE IF UNSPECIFIED IN INPUT
C
IF(TEMP.EQ.0.) TEMP=300.
C
C SPECIFICATION OF THE ERROR BOUND IF UNSPECIFIED IN INPUT
C
IF(ERR.EQ.0.) ERR=1.E-5
C
C OUTPUT OF THE INITIAL CONDITIONS
C
call do_ff(iout)
WRITE(IOUT,5)(ITITLE(I),I=1,7),(REACT(I),I=1,N)
97 WRITE(IOUT,6)(C(I),I=1,N)
C
C COMPUTE THE NET RATES OF REACTION
C
CALL RATES (C,N)
C
C OUTPUT OF THE TEMPERATURE AND ERROR BOUND
C
WRITE(IOUT,8)TEMP,ERR
C
C OUTPUT OF THE DILUTION FACTOR
C
IF(DILUT.NE.0.) WRITE(IOUT,10)DILUT
C
C SET LIMITS FOR TIMED OUTPUTS
C
IF(TSTEP.NE.0.) YMAX(1)=TSTEP
IF(TSTEP.NE.0.) PC=TSTEP
IF(TSTEP.EQ.0.) YMAX(1)=1.E+20
IF(TSTEP.EQ.0.) PC=1.E+20
IF(TFACT.NE.0.) YMAX(2)=TFACT
IF(TFACT.EQ.0.) YMAX(2)=1.
C
C RATE CONSTANTS
C
WRITE(IOUT,9)(R(I),I=1,NR)
GUESS=1.E-16
T=START
IF(NPLOT.EQ.IBLANK) YMAX(3)=0.
IF(NPLOT.NE.IBLANK) YMAX(3)=1.
C
C INITIALIZE PARAMETERS
C
INDEX=1
CALL YFIX (NS,TPRNT,C,NPRNT,INDEX)
C
C USE THE GEAR ROUTINE
C
CALL DRIVES (NS,START,STOPP,C,GUESS,ERR,21,KFLAG)
GO TO 90
C
C CONTINUATION OF DATA
C
98 READ(IN,12)START,STOPP,TEMP,ERR
IF(TEMP.EQ.0.) TEMP=300.
IF(ERR.EQ.0.) ERR=1.E-5
DO 99 I=1,NS
99 REACT(I)=SPECIS(I)
C
C INITIAL CONDITIONS
C
call do_ff(iout)
WRITE(IOUT,5)(ITITLE(I),I=1,7),(REACT(I),I=1,NS)
N=NS
GO TO 97
END
C
C --------------------------------------------------------------------
C
C Output a form-feed character
C
C --------------------------------------------------------------------
C
subroutine do_ff( lun )
character*1 FF
parameter (FF=char(12))
write(lun,'(A)') FF
return
end
C
SUBROUTINE YFIX (N0,TLAST,C,NQ,INDEX)
C
C
C THIS IS THE NORMAL OUTPUT ROUTINE
C
C***********************************************************************
C
COMMON /DATA/ NR, KR( 200,7), A( 200), S( 200), ITITLE(7),
C TEMP, ERR,
C START, STOPP, PC, NP,
C SIG(60),IP(60), ITYPE( 200), R( 200),
C BK, SG, DILUT
COMMON /NAMES/ SPECIS(60), REACT(60), NS
COMMON /STCOM1/N, T, H, HMIN,
C HMAX, EPS1, MF1, KFLAG1,
C JSTART
COMMON /STCOM2/ YMAX(60)
COMMON/INOUT/IN,IOUT,ITAPE,LUNRXR
DIMENSION C(60,6),RT( 200),PA( 60)
2 FORMAT(/,20X,26HSPOT CHECK AT TOTAL TIME =,1PE9.2)
3 FORMAT(/,10H NET RATES,2X,1P 5E12.3 ,/,(12X,1P 5E12.3))
4 FORMAT(//,2X,22HTHE REACTION RATES ARE,/,(12X,1P 5E12.3))
5 FORMAT(/,4X,5HTIME ,4X, 5(4X,A4,4X),/,2X,8HINTERVAL,3X,
C 5(4X,A4,4X),/,(13X, 5(4X,A4,4X)))
C 14 FORMAT(1H1)
16 FORMAT(/,4X,5HTIME ,4X, 5(4X,A4,4X),/)
21 FORMAT(1P 6E12.3/, 6E12.3/,(12X, 5E12.3))
2500 FORMAT(/(1P6E12.3))
C
C ENTER HERE FOR INITIAL CALL
C
GO TO (50,55),INDEX
50 INDEX=2
PA(60)=1.
C
C TIME INTERVALS
C
C INITIALIZE PARAMETERS
C
NT=1
NC=0
ZM=C(NS,1)*1.E-20
ZN=ZM*ERR
NHS=NQ
MS=NQ
TOLD=T
TPRNT=TLAST
TSTEP=PC
ZM=ZM*1.E-14
TFACT=YMAX(2)
IF(TFACT.GT.1.) TSTEP=0.
MT=19-((NS-1)/10)*3-(NR-1)/30
call do_ff(iout)
NN=NS+1-MIN0(NS/10,1)
IF(NN.GE.11) WRITE(IOUT,5)(SPECIS(I),I=1,NN)
IF(NN.LE.10) WRITE(IOUT,16)(SPECIS(I),I=1,NN)
IF(YMAX(3).EQ.0.) NPLOT=0
IF(YMAX(3).NE.0.) NPLOT=1
NTP=0
IF(NPLOT.EQ.0) GO TO 7
IF(T.EQ.START) GO TO 7
GO TO 42
C
C REGULAR ENTRY
C
55 IF(T.EQ.TOLD) RETURN
42 IF(T.GE.TLAST) GO TO 10
TOLD=T
DO 6 I=1,N
C
C CHECK FOR UNREASONABLY LOW CONCENTRATIONS
C
IF(C(I,1).GT.ZN) GO TO 6
C
C INITIATE CALCULATION
C
IF(T.LT.START+EPS1) GO TO 25
C
C IF NEGATIVE CONCENTRATION OCCURS ELIMINATE SPECIES
C
IF(C(I,1).LT.0.) CALL CLEAN (N0,C,I)
C
C CONCENTRATION IS INCREASING
C
IF(C(I,2).GT.0.) GO TO 6
C
C CHECK FOR MINIMUM CONCENTRATION
C
IF(-C(I,2).GT.C(I,1)*ZN.OR.C(I,1).LT.ZM) CALL CLEAN (N0,C,I)
GO TO 6
C
C NEGATIVE CONCENTRATIONS ARE NOT ALLOWED
C
25 IF(C(I,1).LT.0.) C(I,1)=0.
6 CONTINUE
IF(T.GE.TPRNT) GO TO 11
NHS=NHS+1
C
C MS REGULATES THE NUMBER OF PRINT ITERATIONS ALLOWED BEFORE SPOT
C CHECKING THE CONCENTRATIONS WITH THE RATE INFORMATION
C
IF(NHS.GE.MS.OR.T.GE.TLAST) GO TO 7
RETURN
7 NHS=0
NTP=NTP+1
IF(NS.GE.6) WRITE(IOUT,21)T,(C(I,1),I=1,5),H,
C (C(I,1),I=6,NS)
IF(NS.LE.5) WRITE(IOUT,21)T,(C(I,1),I=1,NS),H
C
C MT REGULATES THE NUMBER OF PRINT ITERATIONS PER PAGE
C
IF(NTP.GE.MT) GO TO 8
IF(T.GE.STOPP) GO TO 8
RETURN
8 NTP=0
WRITE(IOUT,2)T
C
C CALCULATE THE NET RATES OF REACTION
C
CALL DIFFUN (N,C,RT)
RT(NS)=0.
DO 85 I=1,N
85 RT(NS)=RT(NS)+RT(I)
WRITE(IOUT,3)(RT(I),I=1,NS)
C
C COMPUTE THE ACTUAL RATES OF REACTION
C
DO 9 I=1,NR
J=KR(I,1)
K=KR(I,2)
L=KR(I,3)
IF(J.EQ.0) RT(I)=0.
IF(J.EQ.0) GO TO 9
IF(J.EQ.99) J=0
IF(J.EQ.0) TJ=1.
IF(J.NE.0) TJ=C(J,1)
IF(K.EQ.0) TK=1.
IF(K.NE.0) TK=C(K,1)
IF(L.EQ.0) TL=1.
IF(L.NE.0) TL=C(L,1)
RT(I)=TJ*TK*TL*R(I)
9 CONTINUE
WRITE(IOUT,4)(RT(I),I=1,NR)
IF(T.GE.TLAST) RETURN
call do_ff(iout)
NN=NS+1-MIN0(NS/10,1)
IF(NN.GE.11) WRITE(IOUT,5)(SPECIS(I),I=1,NN)
IF(NN.LE.10) WRITE(IOUT,16)(SPECIS(I),I=1,NN)
RETURN
C
C INTERPOLATE TO EXACT TIME DESIRED FOR PRINTING OUTPUT
C AT THE LAST TIME STEP
C
10 NHS=MS
NTP=0
TPRNT=TLAST
C
C INTERPOLATE FOR THE EXACT TIME DESIRED TO PRINT
C
11 IF(TPRNT.LT.START) GO TO 15
RR=(TPRNT-T)/H
DO 12 I=1,N
PA(I)=C(I,1)
RH=1.
DO 12 J=1,NQ
RH=RH*RR
12 PA(I)=PA(I)+RH*C(I,J+1)
IF(TPRNT.NE.TLAST) GO TO 95
C
C CALCULATE THE NET RATES OF REACTION
C
95 CALL DIFFUN (N,PA,RT)
IF(NS.LE. 5) WRITE(IOUT,21)TPRNT,(PA(I),I=1,NS),H
IF(NS.GE.6) WRITE(IOUT,21)TPRNT,(PA(I),I=1, 5),H,
1 (PA(I),I=6,NS)
WRITE(ITAPE,2500)TPRNT,(PA(I),I=1,NS)
RT(NS)=0.
DO 80 I=1,N
80 RT(NS)=RT(NS)+RT(I)
WRITE(IOUT,3) (RT(I),I=1,NS)
C
C COMPUTE THE ACTUAL RATES OF REACTION
C
DO 13 I=1,NR
J=KR(I,1)
K=KR(I,2)
L=KR(I,3)
IF(J.EQ.0) RT(I)=0.
IF(J.EQ.0) GO TO 13
IF(J.EQ.99) J=60
IF(K.EQ.0) K=60
IF(L.EQ.0) L=60
RT(I)=PA(J)*PA(K)*PA(L)*R(I)
13 CONTINUE
WRITE(IOUT,4)(RT(I),I=1,NR)
IF(LUNRXR.GT.0)WRITE(LUNRXR,2500)TPRNT,(RT(I),I=1,NR)
NHS=0
IF(T.GE.TLAST) RETURN
NTP=NTP+6
IF(NTP.GT.MT) call do_ff(iout)
IF(NTP.GT.MT) NTP=0
NN=NS+1-MIN0(NS/10,1)
IF(NN.GE. 6) WRITE(IOUT,5)(SPECIS(I),I=1,NN)
IF(NN.LE. 5) WRITE(IOUT,16)(SPECIS(I),I=1,NN)
15 TPRNT=TPRNT*TFACT + TSTEP
RETURN
END
C
SUBROUTINE STIFF (Y,N0)
C
C***********************************************************************
C
C HINDMARSH AT LLL
C
C***********************************************************************
C
COMMON /STCOM1/ N, T, H, HMIN, HMAX,
C EPS, MF, KFLAG, JSTART
COMMON /STCOM2/ YMAX(60)
COMMON /STCOM3/ ERROR(60)
COMMON /STCOM4/ PW(1640),IW(60)
COMMON /STCOM5/ FSAVE(60) ,FSTAR(60)
DIMENSION Y(60,6), EL(6), TQ(4)
DATA EPSOLD/0./, D1/0./, BG/1.E20/
KFLAG=0
TOLD=T
IF(JSTART.GT.0) GO TO 200
NQ=1
IRET=1
L=2
RMAX=1.E4
CRATE=1.
OLDL0=1.
RC=0.
CALL DIFFUN (N,Y,FSAVE)
DO 110 I=1,N
Y(I,2)=FSAVE(I)*H
A=AMAX1(1.,Y(I,1))
110 YMAX(I)=A**(-2)
IDOUB=L+1
HOLD=H
NSQ=N0*N0
130 CALL COSET (NQ,EL,TQ)
RC=RC*EL(1)/OLDL0
OLDL0=EL(1)
140 EDN=(TQ(1)*EPS)**2
E=(TQ(2)*EPS)**2
EUP=(TQ(3)*EPS)**2
BND=(TQ(4)*EPS)**2
GO TO (160,170,200),IRET
150 IF(EPS.EQ.EPSOLD) GO TO 160
IRET=1
GO TO 140
160 EPSOLD=EPS
IF(H.EQ.HOLD) GO TO 200
RH=H/HOLD
H=HOLD
170 RH=AMAX1(RH,HMIN/ABS(H))
RH=AMIN1(RH,HMAX/ABS(H),RMAX)
R1=1.
DO 180 J=2,L
R1=R1*RH
DO 180 I=1,N
180 Y(I,J)=Y(I,J)*R1
H=H*RH
RC=RC*RH
IDOUB=L+1
IF(T.NE.TOLD) GO TO 690
200 IF(ABS(RC-1.).GT.0.3) IWEVAL=1
T=T+H
DO 210 J1=1,NQ
DO 210 J2=J1,NQ
J=NQ-J2+J1
DO 210 I=1,N
210 Y(I,J)=Y(I,J)+Y(I,J+1)
220 DO 230 I=1,N
230 ERROR(I)=0.
M=0
CALL DIFFUN (N,Y,FSTAR)
IF(IWEVAL.LE.0) GO TO 360
CALL PEDERV (N,Y,PW,N0)
R=-EL(1)*H
DO 250 J=1,N
K=(J-1)*N0
DO 249 I=1,N
249 PW(I+K)=PW(I+K)*R
250 PW(J+K)=PW(J+K)+1.
IWEVAL=0
RC=1.
CALL DECOMP (N0,N,PW,IW,FSAVE,IER)
IF(IER.NE.0) GO TO 420
360 DO 370 I=1,N
IF(M.EQ.0) GO TO 365
IF(-H*FSTAR(I)*10. .GT.Y(I,1)) GO TO 420
365 FSTAR(I)=FSTAR(I)*H-Y(I,2)-ERROR(I)
370 CONTINUE
CALL SOLVE (N0,N,PW,FSTAR,FSAVE,IW)
D=0.
DO 390 I=1,N
ERROR(I)=ERROR(I)+FSAVE(I)
D=YMAX(I)*FSAVE(I)**2+D
390 FSAVE(I)=Y(I,1)+EL(1)*ERROR(I)
IF(M.NE.0) CRATE=AMAX1(.9*CRATE,D/D1)
IF((D*AMIN1(1.,2.*CRATE)).LE.BND) GO TO 450
D1=D
M=M+1
IF(M.EQ.3) GO TO 420
CALL DIFFUN (N,FSAVE,FSTAR)
GO TO 360
420 IF(IWEVAL.EQ.-1) GO TO 440
T=TOLD
RMAX=2.
DO 430 J1=1,NQ
DO 430 J2=J1,NQ
J=NQ-J2+J1
DO 430 I=1,N
430 Y(I,J)=Y(I,J)-Y(I,J+1)
IF(ABS(H).LE.(HMIN*1.00001)) GO TO 680
RH=.25
GO TO 170
440 IWEVAL=1
GO TO 220
450 D=0.
DO 460 I=1,N
460 D=D+ERROR(I)*ERROR(I)*YMAX(I)
IWEVAL=-1
IF(D.GT.E) GO TO 500
KFLAG=0
DO 470 J=1,L
DO 470 I=1,N
470 Y(I,J)=Y(I,J)+EL(J)*ERROR(I)
DO 480 I=1,N
YM=Y(I,1)*Y(I,1)
IF(YM.LE.1.E-20) GO TO 480
YT=YM*YMAX(I)
IF(YT.GT.1.) YMAX(I)=1./YM
IF(YT.LT.0.0001) YMAX(I)=1./YM
480 CONTINUE
IF(IDOUB.EQ.1) GO TO 520
IDOUB=IDOUB-1
IF(IDOUB.GT.1.OR.NQ.EQ.5) GO TO 700
DO 490 I=1,N
490 Y(I,6)=ERROR(I)
GO TO 700
500 KFLAG=KFLAG-1
T=TOLD
DO 510 J1=1,NQ
DO 510 J2=J1,NQ
J=NQ-J2+J1
DO 510 I=1,N
510 Y(I,J)=Y(I,J)-Y(I,J+1)
RMAX=2.
IF(ABS(H).LE.(HMIN*1.00001)) GO TO 660
IF(KFLAG.LE.-3) GO TO 640
PR3=BG
GO TO 540
520 PR3=BG
IF(NQ.EQ.5) GO TO 540
D1=0.
DO 530 I=1,N
530 D1=D1+YMAX(I)*(ERROR(I)-Y(I,6))**2
PR3=((D1/EUP)**(.5/FLOAT(L+1)))*1.4+1.4E-6
540 PR2=((D/E)**(.5/FLOAT(L)))*1.2+1.2E-6
PR1=BG
IF(NQ.EQ.1) GO TO 560
D=0.
DO 550 I=1,N
550 D=D+YMAX(I)*Y(I,L)*Y(I,L)
PR1=((D/EDN)**(.5/FLOAT(NQ))+1.E-6)*1.3
560 IF(PR2.LE.PR3) GO TO 570
IF(PR3.LT.PR1) GO TO 590
570 IF(PR2.GT.PR1) GO TO 580
NEWQ=NQ
RH=1./PR2
GO TO 620
580 NEWQ=NQ-1
RH=1./PR1
GO TO 620
590 NEWQ=L
RH=1./PR3
IF(RH.LT.1.1) GO TO 610
DO 600 I=1,N
600 Y(I,NEWQ+1)=ERROR(I)*EL(L)/FLOAT(L)
GO TO 630
610 IDOUB=10
GO TO 700
620 IF((KFLAG.EQ.0).AND.(RH.LT.1.1)) GO TO 610
IF(NEWQ.EQ.NQ) GO TO 170
630 NQ=NEWQ
L=NQ+1
IRET=2
GO TO 130
640 IF(KFLAG.EQ.-9) GO TO 670
RH=10.**KFLAG
RH=AMAX1(HMIN/ABS(H),RH)
H=H*RH
CALL DIFFUN (N,Y,FSAVE)
DO 650 I=1,N
650 Y(I,2)=H*FSAVE(I)
IWEVAL=1
IDOUB=10
IF(NQ.EQ.1) GO TO 200
NQ=1
L=2
IRET=3
GO TO 130
660 KFLAG=-1
GO TO 700
670 KFLAG=-2
GO TO 700
680 KFLAG=-3
GO TO 700
690 RMAX=10.
700 HOLD=H
JSTART=NQ
RETURN
END
C
SUBROUTINE DRIVES (N0,T0,TLAST,Y,H0,EPS,MF,KFLAG)
C
C***********************************************************************
C
C HINDMARSH AT LLL
C
C***********************************************************************
COMMON /STCOM1/ N, T, H, HMIN, HMAX, EPS1,
C MF1, KFLAG1, JSTART
COMMON /STCOM2/ YMAX(60)
COMMON /INOUT/ IN, IOUT, ITAPE
DIMENSION Y(60,6)
100 FORMAT(//44H PROBLEM APPEARS UNSOLVABLE WITH GIVEN INPUT//)
200 FORMAT(//30H KFLAG = -2 FROM STIFF AT T = ,E16.8,5H H =,E16.8/
1 52H THE REQUESTED ERROR IS SMALLER THAN CAN BE HANDLED//)
300 FORMAT(//30H KFLAG = -3 FROM STIFF AT T = ,E16.8/
1 45H CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED/)
C
C N0 IS FOR M GAS
C
N=N0-1
T=T0
H=H0
EPS1=EPS
MF1=MF
KHFLAG=0
JSTART=0
HMAX=AMIN1(YMAX(1),ABS(T0-TLAST)*.1)
HMIN=AMIN1(ABS(H0),.1*HMAX)
10 CALL STIFF (Y,N0)
NQ=JSTART
KGO=1-KFLAG1
GO TO (20,30,50,60),KGO
20 INDEX=2
CALL YFIX (N0,TLAST,Y,NQ,INDEX)
IF((T-TLAST)*H.LT.0.) GO TO 10
GO TO 70
30 IF(KHFLAG.EQ.10) GO TO 40
KHFLAG=KHFLAG+1
HMIN=HMIN*.1
H=H*.1
JSTART=-1
GO TO 10
40 WRITE(IOUT,100)
GO TO 70
50 WRITE(IOUT,200)T,H
GO TO 70
60 WRITE(IOUT,300)T
GO TO 30
70 KFLAG=KFLAG1
T0=T
H0=H
RETURN
END
C
SUBROUTINE COSET (NQ,EL,TQ)
C
C***********************************************************************
C
C HINDMARSH AT LLL
C
C***********************************************************************
C
DIMENSION PERTST(5,3), EL(6), TQ(4)
DATA PERTST /1.,1.,.5,.1667,
*.04167,2.,4.5,7.333,10.42,13.7,3.,6.,9.167,12.5,1./
EL(2)=1.
GO TO (201,202,203,204,205),NQ
201 EL(1)=1.
GO TO 900
202 EL(1)=2./3.
EL(3)=1./3.
GO TO 900
203 EL(1)=6./11.
EL(3)=6./11.
EL(4)=1./11.
GO TO 900
204 EL(1)=.48
EL(3)=.7
EL(4)=.2
EL(5)=.02
GO TO 900
205 EL(1)=60./137.
EL(3)=225./274.
EL(4)=85/274.
EL(5)=15./274.
EL(6)=1./274.
900 DO 910 K=1,3
910 TQ(K)=PERTST(NQ,K)
TQ(4)=.5*TQ(2)/FLOAT(NQ+2)
RETURN
END
C
SUBROUTINE DECOMP (ND,N,U,IS,SC,IER)
C
C***********************************************************************
C
C HINDMARSH AT LLL
C
C***********************************************************************
C
DIMENSION U(ND,N), IS(N), SC(N)
IER=0
DO 1 I=1,N
IS(I)=I
B=0.
DO 2 J=1,N
2 B=AMAX1(B,ABS(U(I,J)))
IF(B.EQ.0.) GO TO 3
1 SC(I)=1./B
IF(N.EQ.1) RETURN
NM1=N-1
DO 5 K=1,NM1
B=0.
DO 6 I=K,N
IP=IS(I)
S=ABS(U(IP,K))*SC(IP)
IF(S.LE.B) GO TO 6
B=S
KP=I
6 CONTINUE
IF(B.EQ.0.) GO TO 4
IF(KP.EQ.K) GO TO 7
J=IS(K)
IS(K)=IS(KP)
IS(KP)=J
7 KP=IS(K)
KP1=K+1
U(KP,K)=1./U(KP,K)
S=U(KP,K)
DO 5 I=KP1,N
IP=IS(I)
U(IP,K)=U(IP,K)*S
B=U(IP,K)
DO 5 J=KP1,N
5 U(IP,J)=U(IP,J)-B*U(KP,J)
IP=IS(N)
IF(U(IP,N).EQ.0.) GO TO 4
U(IP,N)=1./U(IP,N)
RETURN
3 IER=1
RETURN
4 IER=2
RETURN
END
C
SUBROUTINE SOLVE (M,N,RU,B,X,IS)
C
C***********************************************************************
C
C HINDMARSH AT LLL
C
C***********************************************************************
C
DIMENSION RU(M,N), B(N), X(N), IS(N)
IF(N.EQ.1) GO TO 5
NP1=N+1
IP=IS(1)
X(1)=B(IP)
DO 2 I=2,N
IP=IS(I)
K=I-1
S=0.
DO 1 J=1,K
1 S=S+RU(IP,J)*X(J)
2 X(I)=B(IP)-S
X(N)=X(N)*RU(IP,N)
DO 4 K=2,N
I=NP1-K
IP=IS(I)
IP1=I+1
S=0.
DO 3 J=IP1,N
3 S=S+RU(IP,J)*X(J)
4 X(I)=(X(I)-S)*RU(IP,I)
RETURN
5 X(1)=B(1)/RU(1,1)
RETURN
END
C
SUBROUTINE CLEAN (N,X,L)
C
C***********************************************************************
C
C THIS SUBROUTINE ELIMINATES A SPECIES FROM THE SYSTEM
C
C***********************************************************************
C
DIMENSION X(60,6)
COMMON /STCOM4/ PW(1640),IW(60)
COMMON /DATA/ NR, KR( 200,7), A( 200), S( 200), ITITLE(7),
C TEMP, ERR,
C START, STOPP, PC, NP,
C SIG(60),IP(60), ITYPE( 200),
C R( 200),BK, SG, DILUT
DO 1 I=1,N
C
C FIND THE SPECIES NUMBER IN THE LU DECOMPOSED VECTOR
C
IF(L.EQ.IW(I)) GO TO 2
1 CONTINUE
I=N
C
C CLEAR THE APPROPRIATE VALUES OF THE LU DECOMPOSED JACOBIAN
C
2 DO 3 J=1,N
K=(J-1)*N+I
3 PW(K)=0.
DO 4 I=1,NR
DO 4 J=1,7
IF(KR(I,J).NE.L) GO TO 4
IF(J.GT.3) GO TO 5
C
C ZERO ALL RATE CONSTANTS WITH THE SPECIES AS THE REACTANT
C
R(I)=0.
GO TO 4
C
C CHANGE THE SIGN IF THE SPECIES IS A PRODUCT
C
5 KR(I,J)=-L
4 CONTINUE
C
C ZERO THE TAYLOR SERIES OF THE SPECIES
C
DO 6 I=1,6
6 X(L,I)=0.
RETURN
END
C
SUBROUTINE RATES (C,N)
C
C***********************************************************************
C
C THIS SUBROUTINE SETS THE INITIAL CONCENTRATIONS AND CALCULATES THE
C RATE CONSTANTS
C
C***********************************************************************
C
REAL IBLANK,MBLANK,JINTER
COMMON /ALPHA/ IBLANK,MBLANK,JINTER
COMMON /DATA/ NR, KR( 200,7), A( 200), S( 200), ITITLE(7),
C TEMP, ERR,
C START,STOPP, PC, NP,
C SIG(60),IP(60), ITYPE( 200),
C R( 200), BK, SG,
C DILUT
COMMON /NAMES/ SPECIS(60), REACT(60), NS
DIMENSION C(60)
DO 1 I=1,60
1 SIG(I)=0.
DO 3 I=1,N
DO 2 J=1,NS
IF(SPECIS(J).EQ.REACT(I)) SIG(J)=C(I)
IF(SPECIS(J).EQ.REACT(I)) GO TO 3
2 CONTINUE
SIG(NS+1)=SIG(NS+1)+C(I)
3 CONTINUE
N=NS
M=N-1
C(N)=0.
DO 4 I=1,M
C(N)=C(N)+SIG(I)
4 C(I)=SIG(I)
BK=0.
C
C IF SIG(N) DOES NOT EQUAL ZERO, IT IMPLIES THAT THAT THE
C CONCENTRATION OF
C SPECIES N HAS BEEN READ. OTHERWISE M IS THE SUM OF THE INITIAL
C CONCENTRATIONS.
C
C(N)=C(N)+SIG(NS+1)
IF(SIG(N).NE.0.) BK=SIG(N)-C(N)
IF(SIG(N).NE.0.) C(N)=SIG(N)
NP=0
DO 6 I=1,NR
IF(KR(I,1).EQ.0) GO TO 6
ITYPE(I)=2
C
C CHECK FOR 99 CODE AND SUBSTITUTE M WHICH IS THE LAST SPECIES.
C
IF(KR(I,1).EQ.99.AND.KR(I,3).EQ.99) KR(I,1)=N
IF(KR(I,1).EQ. N.AND.KR(I,3).EQ.99) KR(I,3)=0
IF(KR(I,2).EQ.0.AND.KR(I,3).EQ.99) KR(I,2)=N
IF(KR(I,2).EQ.N.AND.KR(I,3).EQ.99) KR(I,3)=0
IF(KR(I,2).EQ.0) ITYPE(I)=1
IF(KR(I,3).NE.0) ITYPE(I)=3
IF(KR(I,3).EQ.99) KR(I,3)=N
C
C RESTORE ANY PREVIOUS CALLS TO CLEAN
C
DO 5 J=4,7
IF(KR(I,J).LT.0) KR(I,J)=-KR(I,J)
5 CONTINUE
C
C CALCULATE THE RATE CONSTANTS
C
R(I)=A(I)*EXP(-S(I)/TEMP)
6 CONTINUE
IF(NS.LE.9) SPECIS(NS+1)=JINTER
RETURN
END
C
SUBROUTINE DIFFUN (N,X,XT)
C
C***********************************************************************
C
C THIS SUBROUTINE CALCULATES THE DERIVATIVE VECTOR OF THE ODE;S
C
C***********************************************************************
C
COMMON /DATA/ NR, KR( 200,7), A( 200), S( 200), ITITLE(7),
C TEMP, ERR,
C START, STOPP, PC, NP,
C SIG(60),IP(60), ITYPE( 200), R( 200),
C BK, SG, DILUT
DIMENSION XT(N), X(N)
X(60)=1.
XT(N+1)=0.
P=BK
DO 1 I=1,N
XT(I)=-DILUT*X(I)
1 P=P+X(I)
X(N+1)=P
DO 2 IR=1,NR
I=KR(IR,1)
IF(I.EQ.0) GO TO 2
C
C CHECK FOR A ZEROTH ORDER REACTION
C
IF(I.EQ.99) I=60
J=ITYPE(IR)
GO TO (11,12,13),J
11 RT=R(IR)*X(I)
GO TO 3
12 J=KR(IR,2)
RT=R(IR)*X(I)*X(J)
XT(J)=XT(J)-RT
GO TO 3
13 J=KR(IR,2)
K=KR(IR,3)
RT=R(IR)*X(I)*X(J)*X(K)
XT(J)=XT(J)-RT
XT(K)=XT(K)-RT
3 IF(I.NE.91) XT(I)=XT(I)-RT
DO 4 K=4,7
I=KR(IR,K)
IF(I.EQ.0) GO TO 2
C
C I WILL BE NEGATIVE IF CLEAN HAS BEEN CALLED
C
IF(I.LT.0) GO TO 4
XT(I)=XT(I)+RT
4 CONTINUE
2 CONTINUE
RETURN
END
C
SUBROUTINE PEDERV (N,X,A,N0)
C
C***********************************************************************
C
C THIS SUBROUTINE GENERATES THE JACOBIAN OF THE ODE;S
C
C***********************************************************************
C
COMMON /DATA/ NR, KR( 200,7), B( 200), S( 200), ITITLE(7),
C TEMP, ERR,
C START, STOPP, PC, NP,
C SIG(60),IP(60), ITYPE( 200), R( 200),
C BK, SG, DILUT
DIMENSION A(N0,N),X(N), RT(3)
C
C INITIALIZE VECTOR
C
X(60)=1.
DO 1 J=1,N0
DO 1 I=1,N0
1 A(I,J)=0.
C
C CHECK FOR NO REACTION OR ZEROTH ORDER REACTION
C
DO 2 IR=1,NR
IF(KR(IR,1).EQ.0.OR.KR(IR,1).EQ.99) GO TO 2
MT=ITYPE(IR)
DO 3 I=1,MT
JINDEX=I+1-I/3*3
KINDEX=I+2-I/2*3
J=KR(IR,JINDEX)
K=KR(IR,KINDEX)
IF(J.EQ.0) J=60
IF(K.EQ.0) K=60
3 RT(I)=R(IR)*X(J)*X(K)
DO 6 K=1,MT
I=KR(IR,K)
DO 4 L=1,MT
J=KR(IR,L)
4 A(I,J)=A(I,J)-RT(L)
DO 5 L=4,7
J=KR(IR,L)
IF(J)5,6,7
7 A(J,I)=A(J,I)+RT(K)
5 CONTINUE
6 CONTINUE
2 CONTINUE
RETURN
END
C
SUBROUTINE MATRX (C,KR)
C
C***********************************************************************
C
C MATRX CREATES A NR X 7 MATRIX OF THE REACTION SCHEME WHERE EACH
C ROW REPRESENTS A REACTION. THE FIRST THREE COLUMNS ARE
C REACTANTS AND THE LAST FOUR ARE THE PRODUCTS. THE ELEMENTS
C CORRESPOND
C TO THE INDIVIDUAL SPECIES AND WILL BE USED AS SUBSCRIPTS.
C
C***********************************************************************
C
REAL IBLANK,MBLANK,JINTER,K,KR
COMMON /ALPHA/ IBLANK,MBLANK,JINTER
COMMON /DATA/ NR, IR( 200,7), A( 200), S( 200), ITITLE(7),
C TEMP, ERR,
C START, STOPP, PC, NP,
C SIG(60),IP(60), ITYPE( 200), R( 200),
C BK, SG, DILUT
COMMON /NAMES/ SPECIS(60), REACT(60), NS
DIMENSION C(60), KR( 200,7)
NOLD=NS+1
DO 10 I=1,NR
C
C SKIP REACTIONS ALREADY PROCESSED
C
IF(IR(I,2).EQ.NOLD.OR.IR(I,3).EQ.NOLD) IR(I,3)=99
IF(IR(I,2).EQ.NOLD) IR(I,2)=0
IF(IR(I,1).EQ.NOLD) IR(I,1)=99
IF(IR(I,1).EQ.NOLD) IR(I,3)=99
IF(IABS(IR(I,1)).NE.100) GO TO 10
C
C IF LESS THAN THREE REACTANTS, FILL FIRST SLOTS.
C
IF(KR(I,1).NE.IBLANK) GO TO 1
IF(KR(I,3).NE.IBLANK) KR(I,1)=KR(I,3)
KR(I,3)=IBLANK
IF(KR(I,1).NE.IBLANK) GO TO 1
KR(I,1)=KR(I,2)
KR(I,2)=IBLANK
1 IF(KR(I,2).NE.IBLANK) GO TO 2
IF(KR(I,3).NE.IBLANK) KR(I,2)=KR(I,3)
KR(I,3)=IBLANK
2 DO 11 J=4,7
C
C GET RID OF M AS A PRODUCT
C
IF(KR(I,J).EQ.MBLANK) KR(I,J)=IBLANK
11 CONTINUE
C
C IF LESS THAN FOUR PRODUCTS, FILL FIRST SLOTS.
C
IF(KR(I,4).NE.IBLANK) GO TO 5
DO 3 J=1,3
INDEX=8-J
IF(KR(I,INDEX).NE.IBLANK) GO TO 4
3 CONTINUE
INDEX=5
4 KR(I,4)=KR(I,INDEX)
KR(I,INDEX)=IBLANK
IF(KR(I,1).NE.IBLANK.OR.KR(I,4).NE.IBLANK) GO TO 5
IR(I,1)=0
GO TO 10
5 CONTINUE
IF(KR(I,5).NE.IBLANK) GO TO 6
IF(KR(I,7).NE.IBLANK) KR(I,5)=KR(I,7)
KR(I,7)=IBLANK
IF(KR(I,5).NE.IBLANK) GO TO 6
KR(I,5)=KR(I,6)
KR(I,6)=IBLANK
6 K=KR(I,6)
IF(K.NE.IBLANK) GO TO 65
KR(I,6)=KR(I,7)
KR(I,7)=K
65 DO 9 J=1,7
K=KR(I,J)
C
C
C PROCESS REACTANTS HERE
C
IF(J.GT.3) GO TO 7
C
C ALL M DEPENDENT REACTIONS ARE TO HAVE A 99 IN THE THIRD SLOT
C
IF(K.NE.MBLANK) GO TO 69
GO TO (66,67,68),J
66 KR(I,1)=KR(I,2)
67 KR(I,2)=KR(I,3)
KR(I,3)=MBLANK
68 IR(I,3)=99
69 K=KR(I,J)
C
C ZERO ORDER REACTIONS HAVE 99 IN FIRST SLOT
C
IF(KR(I,1).EQ.IBLANK) IR(I,1)=99
IF(KR(I,1).EQ.IBLANK) K=99
C
C ALL BLANKS ARE SET EQUAL TO ZERO
C
IF(J.EQ.3.AND.K.EQ.MBLANK) K=99
7 IF(K.EQ.MBLANK.OR.K.EQ.IBLANK) IR(I,J)=0
IF(K.EQ.IBLANK.OR.K.EQ.99) GO TO 9
IF(NS.NE.0) GO TO 12
NS=1
GO TO 13
12 DO 8 L=1,NS
IF(K.NE.SPECIS(L)) GO TO 8
C
C SLOT SET TO SPECIES NUMBER
C
IR(I,J)=L
GO TO 9
8 CONTINUE
C
C IF NO SPECIES ARE FOUND, ADD ONE TO THE LIST
C
IF(SPECIS(NS).NE.MBLANK) NS=NS+1
13 SPECIS(NS)=K
C(NS+1)=C(NS)
C(NS)=0.
IR(I,J)=NS
9 CONTINUE
10 CONTINUE
IF(SPECIS(NS).NE.MBLANK) NS=NS+1
SPECIS(NS)=MBLANK
RETURN
END